Finding Black Hole Apparent Horizons

Black holes are often defined by their event horizons: the surface from which not even light can escape. However, showing that light can never escape means knowing the entire future of the universe. Unsurprisingly, finding the event horizon in an interesting spacetime is usually impractical.

Instead, we often look for apparent horizons. These are surfaces on which, at a "fixed instant" in time, all light is going "inwards". Apparent horizons are local, so can be found by numerical algorithms.

For an in-depth discussion of how to locate apparent and event horizons in numerical simulations, see the Living Review by Thornburg.

Simplified spacetime

In this example we will consider a simple spacetime, as the equations generating the spacetime, and for finding the horizon, simplify greatly. We will assume the spacetime is of Brill-Lindquist type and is also

1. axisymmetric,
2. conformally flat,
3. and time symmetric.

The technical details are less important than the sketch showing the geometry of the problem:

The $z$ axis is the rotation axis. The singularities will sit on the $z$ axis. The apparent horizon(s) will surround the singularities.

With the restrictions above, the spacetime is encoded in a single function $\psi$, which depends on the locations of the singularities. Assume there are $N$ singularities on the $z$ axis at $z^{(i)}$ respectively. Denote the spacetime point at which we want to know the spacetime properties as ${\bf y} = (x, z)$ and the location of the $i^{th}$ singularity as ${\bf y}^{(i)} = (0, z^{(i)})$. Then the single function we need, along with its derivatives, is

\begin{align} \psi & = 1 + \frac{1}{2} \sum_{i=1}^N \frac{m^{(i)}}{\left\| {\bf y} - {\bf y}^{(i)} \right\|_2}, \\ \frac{\partial \psi}{\partial r} & = \frac{1}{2} \sum_{i=1}^N \frac{m^{(i)} \left( z^{(i)} \cos(\theta) - r \right)}{\left\| {\bf y} - {\bf y}^{(i)} \right\|_2^3}, \\ \frac{\partial \psi}{\partial \theta} & = -\frac{1}{2} \sum_{i=1}^N \frac{m^{(i)} r z^{(i)} \sin(\theta)}{\left\| {\bf y} - {\bf y}^{(i)} \right\|_2^3}. \end{align}

Note that we use coordinates $(r, \theta)$, where $(x, z) = (r \cos(\theta), r \sin(\theta))$, to match the symmetries of the spacetime.

Now we can write down the problem that needs solving. The horizon shape function $h(\theta)$, which encodes the radius $r$ of the horizon at the given angle $\theta$, satisfies the boundary value problem

\begin{equation} \frac{d^2 h}{d \theta^2} = 2 h - \frac{\cot(\theta)}{C^2} \frac{d h}{d \theta} + \frac{4 h^2}{\psi C^2} \left( \frac{\partial \psi}{\partial r} - \frac{1}{h^2} \frac{\partial \psi}{\partial \theta} \frac{d h}{d \theta} \right) + \frac{3}{h} \left( \frac{d h}{d \theta} \right)^2, \end{equation}

with boundary conditions

\begin{equation} \frac{d h}{d \theta} \left( \theta = 0 \right) = 0 = \frac{d h}{d \theta} \left( \theta = \pi \right), \end{equation}

should a solution exist. If no solution exists, then there is no horizon. If multiple solutions exist then they do not cross and the AH is the solution with largest $h$ ($r$) for any angle. The auxilliary variable $C$ is defined as

\begin{equation} C = \left( 1 + \left( \frac{1}{h} \frac{d h}{d \theta} \right)^2 \right)^{-1/2}. \end{equation}

We note that $\cot(\theta)$ blows up as $\theta \rightarrow 0$. However, it can be shown that

\begin{equation} \lim_{\theta \rightarrow 0} \frac{\cot(\theta)}{C^2} \frac{d h}{d \theta} = 0, \end{equation}

so that the boundary value problem is well defined for all $\theta$.

Numerical solution

Boundary value problems can often be tackled using the shooting method. This is easily implemented using scipy via scipy.integrate.ode and scipy.optimize.brentq, for example. Here we just display some simple results.

This is a simple problem with a single, common horizon. There is reflection symmetry across the $x$ axis which simplifies the problem further: we can restrict the domain to be $\theta \in [0, \pi/2]$ and impose the boundary condition $dh / d\theta (\theta = \pi / 2) = 0$. This avoids the weak (coordinate) singularity at $\theta = \pi$ (seen in the problematic behaviour of the $\cot(\theta)$ term).

The horizon is not far from spherical in this simple case.

Interactive

Note: the code that follows from here is interactive in some sense, so will not display correctly on the nbviewer website. The interactive 2d plots, where the locations or masses of the singularities may be varied, rely on using version 2.0 of IPython. The 3d plots rely on mayavi. Both must be produced using a local version of the notebook.

We can look at what happens as the singularities move further apart. The horizon becomes steadily more distorted; when the singularities are at $z \approx \pm 0.7$ the horizon develops a "peanut" shape. When the location of the singularities is greater than $z \gtrsim 0.766$ there is no common horizon. The horizon does not "pinch off" - it just no longer exists.

We can also examine the horizon in 3d.

Non-symmetric spacetimes

In the simple examples above, there was a reflection symmetry across the $x$ axis enforced by imposing that every singularity at location $z$ had a singularity of equal mass at location $-z$.

Removing this restriction make the numerical solution via shooting more complex. In the simple case the shooting algorithm guessed the horizon radius $h(0)$ at $\theta = 0$. Together with the boundary condition $dh/d\theta(0) = 0$ this specified an initial value problem; by comparing the result at $\theta = \pi / 2$ we could set up a nonlinear root-finding problem to determine the correct value for $h(0)$. In principle the less restrictive problem only requires us to keep integrating to $\theta = \pi$ and impose the boundary condition there. However, the weakly singular behaviour noted above means that the solution of the initial value problem typically diverges before reaching $\theta = \pi$.

Instead we use shooting to a matching point. We guess the horizon radius at $h(0)$ and at $h(\pi)$. We then integrate the two initial value problems to the matching point at $\theta = \pi / 2$. To be a solution of the boundary value problem, both $h$ and $dh/d\theta$ must match at this point. This gives us a system of nonlinear algebraic equations to determine $h(0)$ and $h(\pi)$.

In addition to being more computationally complex and expensive, it is much more difficult to check that the algorithm fails only when there is no horizon to be found.

It is worth playing with the ranges below and seeing the results. Note that when the algorithm fails, is should be clear from the plot: the surface produced is not sufficiently smooth, for example.

We can again plot the horizon in 3d.